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Abstract 

We use an extension of the diagrammatic rules in random matrix 
theory to evaluate spectral properties of finite and infinite products 
of large complex matrices and large hermitian matrices. The infinite 
product case allows us to define a natural matrix-valued multiplicative 
diffusion process. In both cases of hermitian and complex matrices, 
we observe the emergence of a "topological phase transition", when a 
hole develops in the eigenvalue spectrum, after some critical diffusion 
time Tcrit is reached. In the case of a particular product of two hermi- 
tian ensembles, we observe also an unusual localization-delocalization 
phase transition in the spectrum of the considered ensemble. We verify 
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the analytical formulas obtained in this work by numerical simulation. 
PACS: 05.40.+j; 05.45+b; 05.70.Fh; ll.15.Pg 

Keywords: Non-hermitian random matrix models, Diagrammatic ex- 
pansion, Products of random matrices. 
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1 Introduction 

Applications of random matrix theory (RMT) in physics range from interpre- 
tation of complex spectra of energy levels in atomic and nuclear physics PP, 
studies of disordered systems [3 IS] , chaotic behavior jl] to quantum grav- 
ity OinilZI. Use of progressively developed methods of RMT has also proved 
to be of particular interest in other sciences such as meteorology |HI, image 
processing [HI, population ecology [101 or economy [TTl IT^ . 

In this work, we will be concerned with infinite products of random ma- 
trices (PRM) sampled from general Gaussian ensembles. Products of that 
type appear in various branches of physics and interdisciplinary applications 
[13] where system dynamics is described in terms of evolution operators with 
random coefficients. Perhaps the best known examples are those related 
to thermal properties of disordered magnetic systems, in particular, to the 
localization of electronic wave functions in random potentials [TH 1^. Atten- 
tion has also been focused on applications of PRM to the analysis of chaotic 
dynamical systems [iHl [13 HEl HZ! , where stability of trajectories and their 
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sensitive dependence on initial conditions are measured in terms of charac- 
teristic Lyapunov exponents. Moreover, there are several results exploring 
the use of PRM formalism in applied physics and interdisciplinary research, 
ranging from the studies of stability of large eco- and social- systems [TU] , 
adaptive algorithms and the analysis of system performance under the in- 
fluence of external noises ^H] to image compression ^U] and communication 
via antenna arrays [201 1^ • However, despite the richness of the potential 
applications, the number of analytical results available in the PRM theory is 
still rather limited. 

In this paper we develop effective calculational techniques for studying 
products of random matrices in the large limit, and use these tools to derive 
the properties of the natural matrix valued generalization of the geometric 
(multiplicative) diffusion type process. The scalar versions of these processes, 
leading to log-normal distributions are ubiquitous in various fields [22j and 
we expect the matrix valued generalizations to find numerous applications. 

Let us briefly review the conventional (scalar) multiplicative (geometric) 
random walk in one dimension in the presence of some external drift force. 
The process belongs to the class of (Markovian) Ito diffusion processes, whose 
stochastic variable s undergoes an evolution 

ds 

— = fi{s,t)dt + a{s,t)dWt (1) 
s 

In the above stochastic differential equation (SDE), part of the evolution is 
driven by a variable dWt that represents the Wiener process (integral of the 
Gaussian white- noise) , respecting 

(dWt) = (dW^) = dt (2) 

For simplicity we will limit ourselves to the constant drift fi and constant 
variance a. Finite time evolution of the system could be viewed as a string 
of increments 

(l + /^^ + ^/fxi)(l + /^^ + -\/fx,)... 

■■■(l + /^^ + -\/f^M) 

where M stands for the number of steps in the discretization of the time 
interval. After averaging s{T) over the independent identical distributions 



s(T) = lim 



s{0) (3) 
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(iid) of the Gaussian variables Xj and taking the hmit M oo, we recover^ 
the well known solution for the probability density of s{T): 



p{s,T\so,0) 



exp 



{logjs / So) -fxT+la'Tf 
2a^T 



(4) 



with s restricted to positive values. 

The aim of this paper is to study a similar construction in the space of 
matrices taking the defining equation Q as a guideline for the generaliza- 
tion. We are therefore interested in properties of the matrix- valued evolution 
operator defined as 



Y{T) 



lim 





(5) 



where n is some deterministic "drift" matrix and the stochastic matrices Xi 
belong to identical independent random matrix ensembles. In particular we 
will be studying the eigenvalue distribution of Y{T): 



Pt[z) 



1 

N 



tr5(') {z~Y{T)) 



(6) 



where the average is taken over the stochastic by matrices Xi appearing 
in the definition of Y{T). 

Contrary to several standard multidimensional extensions of Brownian 
walks, we concentrate here on studying how the full spectrum of operator Y 
evolves with time T. 

Let us note some key features of the operator Y. First, since the matrices 
Xi in general do not commute, we are dealing with a 'path ordered product'. 
Second, even if the matrices Xi are hermitian, their product is not, i.e. the 
spectrum in general disperses into the complex plane, showing - as pointed 
out in this paper - some rather unusual feature described as a "topological 
phase transition" . Namely the support of the eigenvalue distribution changes 



^The log-normal distribution s{T) is easy to infer looking at the form of the product 
P|) . Taking the logarithm, expanding and using the central limit theorem we immediately 
see, that the r.h.s. tends to the Gaussian distribution. 
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from a simply-connected two-dimensional island to a two-dimensional island 
with a hole. 

We note that the time evolution of some initial vector |0 > under Y forms 
a very general setup for several multivariate stochastic evolutions of the type 



Analytical results for such processes are however scarce. In a recent study, 
Jackson et al. [221 solved the matrix-valued evolution driven by Gaussian 
noise of the type considered above in the case of 2 by 2 hermitian matrices. 
There are also some known results for the diffusion of the norm of vector |r), 
since this problem can be linked to the Lyapunov exponents of the hermitian 
matrix Y'^Y . We briefly comment on this issue in section 9. 

The main result of this work is the derivation of exact (finite M) and 
asymptotic (M — > oo) spectral formulas for the evolution operators Y , in the 
limit where the dimension N oiY tends to infinity. For simplicity, we set 
to zero the deterministic entries - although one can easily include them 
using our formalism. In this paper we also limit ourselves to ensembles of the 
Gaussian type. Generalizations to non-Gaussian measures (involving finite 
and infinite spectral supports) will be investigated elsewhere. We discuss 
here two cases depending on whether the Xj's are (i) hermitian Gaussian 
matrices (GUE) or (ii) complex (arbitrary) Gaussian matrices (GCE). Similar 
constructions could be used for real and quaternionic matrices. 

The main motivation for this work is to formulate a new and relatively 
simple theoretical framework for studying spectra of the diffusive-like evo- 
lution operators with their further applications to various branches of the 
complex systems analysis in physics, mathematics and in interdisciplinary 
science. 

In particular, the formalism outlined here could be used as a starting 
point for matrix versions of super- or sub-diffusions, generated by e.g. matrix 
Levy analogs of the stable power-tailed distributions. It will be used also for 
the study of spectral statistics of infinite products of unitary operators {e.g. 
Polyakov lines) p4j . 

Most of the derivations presented in this paper will be based on diagram- 
matic methods. In order to make the paper self-contained, in the first two 
sections we briefly present the formalism, first for hermitian RMM (section 2) 
and then in section 3, we recall the extension of diagrammatic methods to 
the non-hermitian case [2Z1 12H] , by using complex Gaussian non- hermitian 




(7) 



6 



random matrix model as an example. In section 4 we move to the main 
part of our paper and introduce our method for the treatment of products of 
random matrices by performing calculations for a product of two matrices. 
Then (section 5) we move toward the products of M complex matrices, and 
we perform the hmit M — > oo. We present the final formulas for the spectral 
density and the equations defining the boundary of the two-dimensional sup- 
port of the eigenvalues. In sections 6 and 7 we consider the case of hermitian 
matrices which, surprisingly, is more difficult than the previous one. We an- 
alyze first (section 6) the product of two matrices. This is interesting on its 
own, since it exhibits quite singular behavior - for evolution times smaller 
than a certain critical value, the a priori complex eigenvalues are "frozen" 
to the real axis. In section 7 we move toward the general case. For prod- 
ucts of matrices involving more than two hermitian matrices, the eigenvalues 
always form two-dimensional islands. We solve the spectral problem in the 
large M limit. In section 8 we show the numerical analysis, confirming our 
analytical predictions and the existence of "topological phase transitions", 
corresponding to structural changes of the spectrum. 

Finally, in section 9 we briefly discuss the problems of stability of such 
products and the relation of our results to the standard Lyapunov exponents. 
Section 10 concludes the paper. 

2 Hermitian Random Ensembles 

A key problem in random matrix theories is to find the distribution of eigen- 
values Aj, in the large N (size of the matrix H) limit, i.e 



where the averaging (...) is done over the ensemble oi N x N random her- 
mitian matrices generated with probability 



The eigenvalues of course lie on the real axis. By introducing the resolvent 
(Green's function) 





P{H) oc e 



NTrV{H) 



(9) 




(10) 
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with z = zIn and by using the standard relation 



P.V.\^i7i6{\) . (11) 



X±ie A 

the spectral function (jS)) can be derived from the discontinuities of the Green's 
function ^ 

^lim(G(A-^e)-G(A + ^6)) = I {Tr 6{X-H)) 



\i=i 



N \ 

E'5(A-A.)) = p(A). (12) 



There are several ways of calculating Green's functions for HRMM ^ [3 
0]. We will follow the diagrammatic approach introduced by A starting 
point of the approach is the expression allowing for the reconstruction of the 
Green's function from all the moments (Trif"), 



z z z z z z 



= ^E^(M") (13) 

The reason why the above procedure works correctly for hermitian matrix 
models is the fact that the Green's function is guaranteed to be holomorphic 
in the whole complex plane except at most on one or more 1-dimensional in- 
tervals. We will use the diagrammatic method to evaluate efficiently the sum 
of the moments on the right hand side. We will now restrict ourselves to the 
well known case of a random hermitian ensemble with Gaussian distribution. 

The first step is to introduce a generating function with a matrix-valued 
source J: 

where we integrate over all A^^ elements of the matrix H. All the moments 
follow directly from Z{J) through the relation 



(M") = ^Tr(-) Z{J)^^^ (15) 
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Figure 1: Large N "Feynman" rules for "quark" and "gluon" propagators. 




Figure 2: Diagrammatic expansion of the Green's function up to the 0{H^) 
terms. 

and are straightforward to calculate, since in the Gaussian case the partition 
function (fT^ reads Z{J) = exp j^TtJ"^ . Accordingly, the lowest nonzero 
expectation value is 

and the next non-vanishing expectation value reads 

(H^mm) = ^ {sim^i + ^t^'M + s^fS'Mt) (17) 

The key idea in the diagrammatic approach is to associate to the expressions 
for the moments, like the above, a graphical representation following from a 
simple set of rules. The power of the approach is that it enables to perform 
a resummation of the whole power series (fTSj) through the identification of 
the structure of the relevant graphs. 

We depict the "Feynman" rules in Fig. ^ similar to the standard large 
diagrammatics for QCD jSEI- The 1/z = l/z6^ in is represented by a 
horizontal straight line. The propagator (fT^ is depicted as a double line. 

The diagrammatic expansion of the Green's function is visualized in 
Fig. 121 where one connects the vertices with the double line propagators 
in all possible ways. Each "propagator" brings a factor of and each 

loop a factor of 6^ = N. From the three terms, corresponding to (jl7|) con- 
tributing to (tr H^) only the first two are presented in Fig. |21 (the third and 
the fourth diagram). The diagram corresponding to the third term in ()17p 
represents a non-planar contribution which is suppressed as l/N"^ and hence 
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Figure 3: Schwinger-Dyson equation for rainbow diagrams. 

vanishes when N ^ oo. In general, only planar graphs survive the large 
limit. 

The resummation of (fT!^ is done by introducing the self-energy S com- 
prising the sum of all one-particle irreducible graphs (rainbow-like). Then 
the Green's function reads 

In the large limit the equation for the self energy S, follows from 
resumming the rainbow-like diagrams of Fig. |21 The resulting equation 
( "Schwinger-Dyson" equation of Fig. ^ encodes pictorially the structure of 
these graphs and reads 

E = ^TrGl = ^G = G. (19) 

Equations (fT^ and (fTIHl give immediately G{z — G) = l which can be solved 
to yield 

1 



G^(z) = -izTV^^) (20) 

Only the G- is a normalizable solution, with the proper asymptotic behavior 
G{z) — >• 1/z in the z — > cxd limit. From the discontinuity (cut), using (O, 
we recover Wigner's semicircle for the distribution of the eigenvalues for 
hermitian random matrices 



p{X) = (21) 
Zn 



3 Non-hermitian Random Ensembles 

The main difficulty in the treatment of non-hermitian random matrix models 
is the fact that now the eigenvalues accumulate in two-dimensional domains 
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in the complex plane and the Green's function is no longer holomorphic. 
Therefore the power series expansion (|13|) no longer captures the full infor- 
mation about the Green's function. In particular the eigenvalue distribution 
is related to the non-analytic (non-holomorphic behaviour of the Green's 
function: 

-d,Giz) = piz) . (22) 

71 

This phenomenon can be easily seen even in the the simplest non-hermitian 
ensemble — the Ginibre-Girko one pUj I31j. with non-hermitian matrices X, 
and measure 

P(X) = e-^^^^' . (23) 

It is easy to verify that all moments vanish (trX") = 0, for > so the 
expansion (jl3p gives the Green's function to be G{z) = 1/z (diagrammati- 
cally this follows from the fact that the propagator < X^^X^ > vanishes and 
hence the self-energy S = 0). The true answer is, however, different. Only 
for \z\ > 1 one has indeed G{z) = 1/z. For \z\ < 1 the Green's function is 
nonholomorphic and equals G{z) = z. 

In spite of the above difficulties a diagrammatic approach can be used 
for recovering the full information including the eigenvalue density p{z). The 
first step is to define a regularized Green's function from which one can ex- 
tract p{z). This has been done by exploiting the analogy to two-dimensional 
electrostatics |31|^E2]- The resulting regularized Green's function is how- 
ever difficult to calculate. The second step is to give an equivalent linearized 
form which can form a basis for a diagrammatic calculation j2Z]. We will 
now briefiy review the above developments. 

Let us define the "electrostatic potential" 

" Trln[(z-7W)(z-7W^) + e2]. (24) 



Then 



N 



lim^-— ^— = lim— (Tr- 



e-o dzdz e^oN\ (Iz-A^P+e 



2^2 , 



^(y:S^'\z-K))^npix,y) (25) 
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represents Gauss law, where z = x + iy. The last equality follows from the 
representation of the complex Dirac delta 



e2 



6^'\z - A.) = hm , ,^ , (26) 



-0 (e2 + \z-\ 



In the spirit of the electrostatic analogy we can define the Green's function 
G{z,z), as an "electric field" 

Then Gauss law reads 

d,G = npix,y) (28) 



in agreement with 

As mentioned above, instead of working ab initio with the object ()27|). 
and in view of applying diagrammatic methods it is much more convenient 
to proceed differently - as we will now discuss. 

Following [27j we define the matrix- valued resolvent through 



1 / (z-X le ' 



with 



A 
B 
C 



(z-Xt)(z-X) + e2 
—ie 

(z-X)(z-Xt) + e2 

—ie 

(z-Xt)(z-X) + e2 



and where we introduced the 'block trace' defined as 
^ ( A B\ _/TrATr5\ 

^'^n^^J..><.rlTrCTrDj^;^ (^1) 
12 



Then, by definition, the upper-right component Qu, is equal to the Green's 
function (ITTfl . 

The block approach has several advantages. First of all it is linear in 
the random matrices X allowing for a simple diagrammatic calculational 
procedure. Let us define 2N by 2N matrices 

Then the generalized Green's function is given formally by the same definition 
as the usual Green's function G, 

What is more important, also in this case the Green's function is completely 
determined by the knowledge of all matrix- valued moments 

Ttb2 ^-^nZ-^n.-.Z-^) . (34) 



This last observation allows for a diagrammatic interpretation. The Feyn- 
man rules are analogous to the hermitian ones, only now one has to keep 
track of the block structure of the matrices, e.g. single straight lines will 
now be associated with a matrix factor Z~^. We will now demonstrate the 
above procedure by solving diagrammatically the complex Gaussian Random 
Matrix Model. 

Example: the Girko-Ginibre ensemble 

The ensemble is defined by the measure 

P{X) oc e-^T'^^^' . (35) 
In this case the double line propagators are 

{X^X^D = {xKx^,) = j^5^A- (36) 

As previously, we introduce the self-energy S, (but which is now matrix- 
valued), in terms of which we get a two by two matrix expression 

g = {Z-t)'\ (37) 
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where the 2 by 2 matrix Z reads 



Z= \ ^ (38) 



le z 



The resummation of the rainbow diagrams for S is more subtle, but follows 
easily from the structure of the propagators (jHBj) . The analogue of (fT^ is 
now: 









r 








1 = 


Wii 


j 



E = = • (39) 



The two by two matrix equations (j37ll39p completely determine the prob- 
lem of finding the eigenvalue distribution for the Girko-Ginibre ensemble. 
Inserting ^ into ^ we get: 



Qii Qii \ ^ 1 / z g^j 

Gil ) - ^li^ii V ^Ti ^ 



(40) 



Note that at this moment we can safely put to zero the regulators e. Looking 
at the off-diagonal equation 

Gn = I I, r (41) 

we see that there are two solutions: one with Qii = 0, and another with 
Qil 7^ 0. The first one leads to a holomorphic Green's function, and a 
straightforward calculation gives 

G{z) = i . (42) 

The second one is nonholomorhic, imposing the condition 

\z\^~b^ = l (43) 

where we denoted QiiQii = b"^, hence 

G{z, z) = z (44) 

which leads, via the Gauss law, to 

p{x,y) = ^^_G{zrz) = -. (45) 
14 



Both solutions match at the boundary 6^ = , which in this case reads zz = 1. 
In such a simple way we recovered the results of Ginibre and Girko for the 
complex non-hermitian ensemble. The eigenvalues are uniformly distributed 
on the unit disk \z\'^ < 1. 

This example illustrates more general properties of the matrix valued gen- 
eralized Green's function. Each component of the matrix carries important 
information about the stochastic properties of the system. There are always 
two solutions for Qn, one holomorphic, another non-holomorphic. The sec- 
ond one leads, via Gauss law, to the eigenvalue distribution. The shape of 
the "coastline" bordering the "sea" of complex eigenvalues is determined by 
the matching conditions for the two solutions, i.e. it is determined by im- 
posing on the non-holomorphic solution for the equation 6^ = 0. For other 
important features of the explained above diagrammatic method and more 
complex examples see ES] • 



4 Product of Two Complex Matrices 

We will now move toward the main goal of this paper, namely the evaluation 
of the spectral properties of the random matrix products of the form (0) - 
with zero drift and normalized variance: 



Y{T) 



hm Ym 

M^oo 



lim 




(46) 



In principle the methods of section 3 could be used to study such random 
matrices, however the nonlinear product structure would make the resulting 
diagrammatic rules exceedingly difficult to control. Fortunately, for the pur- 
poses of calculating the Green's function, and hence the eigenvalue density, 
the problem may be linearized at the cost of increasing the size of the matri- 
ces. In this section we will show how the method works for Y2 - the product 
of two matrices, while in section 5 we will consider the general case including 
the limit M —>■ 00. We therefore have to calculate 



Y, 



1 + 



-X2 ) = A,A2 



(47) 
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where Xi and X2 belong to identical independent Gaussian Complex En- 
sembles (iiGCE=Girko-Ginibre). Since the matrices are non-hermitian we 
have to use the formalism of section 3 with a conjugate 'antiholomorphic' 
copy. /^From now on in the rest of the paper we set the regulator e to zero. 
In order to linearize the product structure we use a trick: we introduce a 
further doubling, thus constructing an auxiliary AN by AN Green's function 



( 



w 








w 











\ 


/ 





A 










A2 





w 


















w 


/ 


V 









/ 


w 


-1 








\ 




-1 


w 
















w 


-1 




V 








-1 


w 





-/2X 







A\ 


X2 









A 


Xi 






-1 

















ANxiN 



(48) 



In the above we first separated the "deterministic" part from the "random 
one" (second line) and in the last equality we introduced the notation similar 
to the hermitian and non-hermitian cases analyzed in the previous sections. 

In analogy to (jHT|) . we define now a block-trace operation tr 54, where we 
trace each by block of the AN by AN matrix Q{w) separately. In such 
a way, we obtain a 4 by 4 auxiliary Green's function 



9{w) 



( 9ii 


9l2 


9ii 


9l2 \ 


921 


922 


921 


922 


9ii 


9l2 


9ii 


9l2 


\ 921 


922 


921 


922 / 



— tr B4G{w) 



(49) 



4x4 



The advantage of this function lies in the fact, that the eigenvalues of the 
product A1A2 coincide with the squares of the eigenvalues of the 2N x 2N 
block matrix 

( A,\ 
[ A2 



(50) 



This is due to the off-diagonal structure of the sub-blocks in (j^UI). 
general discussion for an arbitrary product see section 5. 



For a 
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As before we define a 4 by 4 matrix of self-energies 

g{w) = {yv-t)-^ 



(51) 



where W is the result of block-tracing 'bold' W. 

Similarly as in the Girko-Ginibre case, we can now diagrammatically an- 
alyze the content of the matrix S. We see that only the "double line" prop- 



agators (XiX[) and 



X2X2) are different from zero. 



Therefore from 



all of the 16 elements of matrix S only four are nonzero. Moreover, due to 
identical measures of Xi and X2 variables, we get 



^22 

^22 



ag22 
ag22 



(52) 



where the factor a = t/2 comes from the propagator in the corresponding 
Schwinger-Dyson equation of the type discussed in the previous chapter. 

Filling the matrix S with entries (j52p we arrive at the 4 by 4 matrix 
equation for the elements of the Green's function: 



( gii gi2 gn gu \ 

g2i g22 g2i g22 

gu gi2 gn gn 

V g2i fi'22 g2i fl'22 J 



4x4 



W 





V 
/ 



w 








w 






c^gn 









ctgn 



\ 


-1 

w / 

agii 

agii 

/ 



(53) 



As in the Ginibre case, we solve it for gii and gii. Then, using the 
solution, we calculate 5^11. Coming back to original variables with w"^ = z we 
set G{z,z) = gii{w). The eigenvalue distribution then follows from 



1 w _ 
pix,y) = d^G{z,z) 

IT Z 



(54) 



This is a consequence of the relation between the eigenvalues of the block 
matrix (jMIj) and the eigenvalues of the product A1A2. A general derivation 
is outlined in the following section. 

A little algebra (basically an inversion of the 4 by 4 matrix) shows, that, as 
before, we get two solutions for 6^ = giigii= 5'225'22- (i) the holomorphic one 
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Figure 4: Evolution of the contour of the non-holomorphic domain in the 
eigenvalue spectrum of the product of two Ginibre-Girko ensembles as a 
function of evolution times r = 0.1,0.25,0.5,1,1.5,2,3,4,5,6, rising from 
inside to outside. 



(corresponding to 6^ = 0, i.e. all equal to zero); (ii) the non-holomorphic 
one, given by the equation 

^ {\w\^ - a%'' + l) = {a%'' - 1 + w''){a%'' - I + w^) - a%\w + (55) 

The boundary of the eigenvalue support is given, as before, by the above 
equation with 6^ set to 0. It is therefore a fourth order conchoid-like planar 
curve, which in polar coordinates is given by 

-(1 + r) = + 1 - 2r cos(/) (56) 

On Fig. m we show the shape of the nonholomorphic island for several 
sample "evolution times" r. 

The Green's function for the product follows from 

where 

det = {a%'^ - 1 + w'^){a%'^ - 1 + w^) - a%'^{w + (58) 
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5 Product of Arbitrary Number of Complex 
Matrices 



We consider now the product of an arbitrary number of matrices. To see the 
pattern, let us look briefly at the case of three matrices 

n = (i + /Ix.) (l + f^X,) (l + f^x, 

= AiA^A, (69) 

where Xi, X2, X2 again belong to the Girko-Ginibre ensemble. 

Our approach again follows from a very simple exact relation between the 
eigenvalues of a product of N x N matrices 741^42^43 and the eigenvalues of 
a block matrix 

/ Ai \ 
83=1 A2 \ (60) 

W3 ;3^^3^ 

Indeed if {Xi} are the eigenvalues of A1A2A3 then the eigenvalues of the block 
matrix are {X]^^ , X]^^ ■ e^, Xy^ ■ e^}. This is an exact relation for any N 
and follows from the relation between the resolvents (here the matrices Ai 
are of finite size and fixed i.e. no averaging) 

^'■^ = 3^'^ ^ ^^^^^^^^ = z - MA2A, 

namely 

wGb^{w) = zGa^A2A3{w^ = z) (62) 

This is due to the cyclic structure of the block matrix (jUUj) . Obviously, only 
multiplicities of the cubic powers of ^3 contribute to the trace. 

The relation between the eigenvalues now follows from the location of the 
finite number of poles of both functions. 

We can thus safely calculate the eigenvalue density for the block matrix 
S3 using the diagrammatic methods and then extract the density for the 
product through 

Pa,A2...Am{z, z) = ^d2GA,A2...AM = J^-^P^Aii'^, ^) (63) 
where M = 3 and w*^ = z and pism{'Wj'^) = ^du,GB{w,w). 
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To extract only the moments im^olving powers of the triples we 
construct an auxiliary 6A^ by 6A^ Green's function, triplicating the A1A2A2 
products by rewriting them as cyclic block matrices: 
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(64) 



6Nx6N 



where again we separated the "deterministic" part from the "random one" . 
Introducing now the block-trace operation tr b6, we obtain a 6 by 6 auxiliary 
Green's function g{w) = tr b6Q{w). 

The generalization for arbitrary M is now straightforward. For 



M 



(65) 



we define an auxiliary 2MN by 2MN Green's function of the form 



(66) 



where the blocks 



/ w -1 ... \ 

w -1 ... 

... w -1 

\-i ... w y 



(67) 



MNxMN 
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and 



f Xi 
X2 



X 
















Xm-1 




(68) 



MNxMN 



are themselves NM by NM matrices, i.e. each of the hsted elements in W 
and in X is itself an by matrix, either diagonal, denoted by a bold 
symbol, or a random entry Xi, otherwise an by block of zeroes. 

We take now a block-trace operation tr bm, where we trace each A^ by 
N block of the 2MN by 2MN matrix Q{w) separately. In such a way, we 
obtain an 2M by 2M auxiliary Green's function 
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As before we define a 2M by 2M matrix of self-energies Ejj 

-1 



9[w) 



W 

wt 



-E 



(70) 



where W is a result of block-tracing W. We can now diagrammatically 
analyze the content of the matrix E. As previously, only "double line" prop- 
agators (^XiXj^^ for i = 1, . . . , M are different from zero. Therefore from 

all of the 4M^ elements of the matrix E only 2M are different from zero. 
Due to the symmetries we get 



^MM = «5'll = (^9lM = OLQ 



(71) 



where the factor a — t/M comes from the propagator in the corresponding 
set of Schwinger-Dyson equations as in the previous chapter. 
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Inserting now the matrix S with entries (|71|) we arrive at a 2M by 2M 
matrix equation for the elements of the Green's function. First, we solve it 
for g and g. Second, using the solutions, we calculate gu. Third, using (jUHj) 
we get an explicit equation for the spectral density 

, , 11 WW ^ , 
p(a;, y) = -Y7—dw9u{w, w) (72) 

71 IVl ZZ 

where z = x ^ iy. For arbitrary M, the algebra is rather involved. Luckily, 
both the block and the cyclic structure of the main entries W and S make 
taking an inverse of the 2M by 2M matrix possible. The inverse of a cyclic 
matrix is a cyclic matrix, and its explicit form can be obtained from the 
solution of an associated recurrence relation, e.g. using the transfer matrix 
techniques which we will now proceed to do. 

At this moment, it is useful to introduce the notation 

1 + |wp — = 2|w|coshM 

\w\lw = 6 (73) 

The inversion needed for the Green's function reduces to finding explicit 
forms of the cyclic M by M matrices C, C defined as 

C = [W ■ - a^H]-' C =[W^ - W -a%H]'' (74) 

where 6^ = gg. C is just the transpose, or equivalently, the complex conjugate 
of C. The elements cq, . . . cm-i of the first row of the cyclic matrix C fulfill 
the recurrence relation 





(75) 



with the 'boundary conditions' 






(76) 

The above recurrence can be easily solved using transfer matrix techniques, 
i.e. by diagonalizing the 2 by 2 transfer matrix in (fTSj) and returning at 
the end of the calculations to the original basis. An explicit solution for the 
cyclic string c^, for k = 0, 1, .., .M — 1 reads 

= 1 / At _ A^, \ 

w(A+-A_) I Ai^-1 A*:^-li ^ ^ 
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where A± = Se^'^ denote the eigenvalues of the transfer matrix. 

We may now express our key quantities for arbitrary M in terms of ((TTj): 

guiw) = [W^ ■ C]ii = w ■ Co - Cm_i 
g = ag[C]u = ag-cl 

g = ag[C]ii = ag-co (78) 

where * denotes complex conjugation. As always in the case of the complex 
spectra, we have two solutions: trivial (holomorphic), corresponding to g = 
g = 0, and the nontrivial (non-holomorphic), corresponding to the case 6^ = 

Inserting the nonholomorphic solution into gn, redefining G{z, z) = gii{w) 
with = z, and using yields the eigenvalue density for Ym- Since 
du,z = Mz/w, the result for the final spectral function can be further simpli- 
fied 

p(x,y) = --d,Giz,z) (79) 

7T Z 

Explicit formulas for arbitrary M could be easily constructed e.g. by 
running any symbolic program using our solution. The shape of the bound- 
ary is given by the condition that the nonholomorphic solution meets the 
holomorphic one on the complex plane. 

After solving the problem for arbitrary M, we can now address the orig- 
inal problem, i.e. the solution for a continuous product of matrices, corre- 
sponding to the limit M ^ oo. This means, that we have to perform a careful 
limiting procedure in the above formulas. Introducing U = Mu, expanding 
in M, and using the polar parametrization z = r exp i(f), we get the solution 

^, , 1 1 ill sin (h , . 

G(.,.) = -logr + ----^ (80) 



where U = y log r — r^fo^ with b fulfilling the equation 

cos = cosh t/ — sinh [/ (81) 

The set of equations ()80I81|) completes the solution of the continuous product 
of large complex matrices with Gaussian disorder. We would like to note. 
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that due to the hmiting procedure M — oo differentiation of the Green's 
functions has to be done with care. Taking the large M hmit of (j79j) yields 

p{x,y) = --d,G{z,z) (82) 

TT Z 

As a cross-check of our calculations we verified that the spectral function 
obtained in this way is real. 

Substituting 6^ = in (|HT|). we get the equation for the shape of the sup- 
port of the eigenvalues for arbitrary r. Note that the shape of the boundary 
depends on the variable log reflecting in the spectrum the multiplicative 
(geometric) feature of the stochastic process induced by the multiplication 
of random matrices. Explicit solution for the boundary reads 

- =r2 + l-2rcos0 (83) 

2 \ mr J 

It is rewarding to compare solution ()83p to the solution of the conchoid-like 
boundary in the case of two complex matrices ()56|) . In the limit of infinitely 
many products, the shape of the boundary acquires a new symmetry - (j83|) is 
invariant under inversion operation r — > 1/r. It is visible also from the gen- 
eral form ()8H). since the radial dependence for 6^ = is a function of | logr| 
only. This symmetry is responsible for the appearance of a structural change 
of the spectrum - provided r is sufficiently large, two boundaries, related by 
inversion in the radial variable, appear. We describe this structural change of 
the complex spectrum as a "topological phase transition". In Section 8, we 
study numerically the spectral density predicted by our formulas, in partic- 
ular, the "diffusion" of the shape of the boundary caused by evolution with 
the "time" parameter r. We also confirm numerically the critical behavior 
for certain Tcrit. 



6 Product of Two Hermitian Matrices 

In this section we will consider the spectral distribution of the product 

Y,= (l + yiiJi) (l + ^H,^ (84) 

where Hi,H2 are not arbitrary complex matrices, but they are identical in- 
dependent hermitian Gaussian matrices (iiGUE). Contrary to naive expec- 
tations, the above problem is more subtle then the analogous problem of the 
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two complex Gaussian matrices. First, let us note, that the product of two 
hermitian matrices is not, in general, a hermitian matrix, so the spectrum 
of the strings of hermitian matrices develops complex eigenvalues. 

Second, the case of only two hermitian matrices of the considered type 
turns out to be in some sense singular - for certain values of the evolution 
parameter r, the a priori complex eigenvalues get localized on the real axis. 
This phenomenon reminds a bit the non-hermitian localization in the Hatano- 
Nelson model jHEI- Third, the problem is algebraically more involved since 
there will be more nonzero propagators than in the complex case. 

Since even a product of two hermitian matrices can develop a priori a 
complex spectrum, we have to use the same type of auxiliary Green's func- 
tion as in the complex Gaussian case. The construction parallels exactly the 
scheme for two complex matrices considered previously, with the matrices Xi 
and X] replaced both by Hi. The first new feature appears when we con- 
struct the "Schwinger-Dyson" equations for the 4 by 4 self-energy matrix E. 
Since the matrices are hermitian [Hi = Hj), additional nonzero contractions 
appear in the matrix E 
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(85) 



where again the factor a = r/2 comes from the propagator, h and h are 
the new entries and we have exploited the symmetries of the problem (e.g. 
5^11 = 922 etc). 

The solution of the 4 by 4 matrix equation for the 4 by 4 Green's function 

-1 



g{w) = (W-E) 



(86) 



is more involved, since we have to solve it not only for the "gaps" g and g but 
also for additional gaps h and h. Only then we can calculate gu, to recover 
the wanted eigenvalue distribution in the same way as in the complex case 

(El- 
Introducing (P = [1 + ah){l + ah) and = gg, we obtain two (modulo 
their tilded partners) gap equations 

det 



9 
h 



— (a^r(l + ah) + (1 + a/i)(w^ - (1 + ah)' 



(87) 
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Figure 5: Complex eigenvalues of 10 products of two 200 x 200 hermitian 
matrices before and at the transition {\Jt/2 = 0.4 and 0.5). The x and y 
axis correspond to 3?A and respectively. 



-o.H 

Figure 6: Eigenvalues of 10 products of two 200 x 200 hermitian matrices 
after the transition (y'r/2 = 0.7 and 1.0). The axis are defined as in the 
previous figure. 
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coupled to the needed gu 

= ^[-a^b^w + w{w^ - {1 + ahf) (88) 

where det is the determinant of the matrix W — S and reads explicitly 

det = aV + {w^ - (1 + ahf){w^ - (1 + ah)^) - 2a%\d^ + \w\'') (89) 

Numerical simulations show that the eigenvalue distribution for the product 
of two matrices is concentrated on an interval on the real axis for r < 1/2 and 
only then when the edge of the interval reaches the origin, do the eigenvalues 
start to spread into the complex plane (see figures and IH)) • Let us show 
how to find the critical value r = 1/2 from the above formulas based on the 
above numerical observation. 

We thus have to find the point when the boundary of the two dimensional 
island hits the origin w = Then the equations (|H7|) and (|HI?|) simplify to 

h = , (90) 

1 + ah ^ ^ 



26 



1 = ahh (91) 

where we set = gg = 0. Then it is easy to verify that this set of equations 
is satisfied for a = 1/4 which is equivalent to r = 1/2. 



7 Products of Arbitrary Number of Hermi- 
tian Matrices 

Finally, we consider an arbitrary number of matrices belonging to iiGUE. 
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(92) 



As before, the construction of the auxiliary Green's function parallels the 
complex case, up to the point when we have to construct explicitly the ma- 
trix of self-energies E, corresponding to the set of 4M^ Schwinger-Dyson 
equations. Generalization of the results from the previous sections leads to 
the explicit form of the auxiliary Green's function. Introducing 



W S 



(93) 



with blocks 



and 



w 




—1 — ah 



-1 — ah 

w 






S = —a 





-1 — ah 



9 
g 




V 



w 




\ 







-1 — ah 

w 



(94) 



MxM 



9 

9 / MxM 



(95) 



27 



the consistent set of equations comes as 



9ii = Gi,i 
h = Gi,M 
h = Qm+i,m+2 

9 = Gl,M+l 

9 = Gm+1,1 (96) 

i.e. requires an exphcit inversion of the 2M by 2M matrix. The algebra is 
more comphcated since, hke in the case of two hermitian matrices, we have 
now two types of gap equations (for g and h and their tilded partners). Using 
the solutions, we calculate gu, then, we get an explicit equation for the spec- 
tral density (jUHj) . Finally, we perform the limiting procedure. Technically, 
one has to repeat the transfer matrix technique introduced in section 5 for 
the complex case. The algebra is lengthy, but due to the cyclic nature of 
the blocks of matrix TZ the results form a rather surprisingly simple set of 
equations. Returning to the original variable z = rexp(20) we get 

1 3 il-f 

G{z,z) = — logr H sin-?/'/ sinh [/ (97) 

2r 4 r 

and the gap equations read 

/i = — log ^ ^ ~ " — sin -ip I sinh \J = G — \ (98) 



and 



T 

COS '\\) = cosh U — —— sinh U (99) 



where ip = (p + and U = ^ + r/Af - r^gg. The above set of 

equations completes the solution of the continuous product of matrices with 
Gaussian hermitian disorder. Spectral density follows from (j82p . Note, that 
the dependence on log \ z\ reflects the multiplicative character of the Brown- 
ian motion. Again, the boundary has also an inversion-type symmetry, but 
now it is of the form \z\ — > exp(— r)/|z|. Due to the non-zero gap h, the 
boundary in the hermitian case is given by the pair of coupled transcenden- 
tal equations and can be presented numerically only. The appearance of the 
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aforementioned symmetry is responsible for the "topological phase transi- 
tion" in the spectrum of the product of infinitely many matrices, similar to 
the one observed for the product of infinitely many complex matrices. 

In Section 8, we study numerically the spectral density predicted by our 
formulas, in particular, the "diffusion" of the shape of the boundary caused 
by evolution with the "time" parameter r. We confirm the appearance of a 
"phase transition" for certain Tcrit in the hermitian case as well. 



In this section we present the results of numerical simulations, confirm- 
ing our analytical predictions. We start from the complex case. Figure [7| 
shows the evolution of the boundary for several sample evolution times 
r = 0.5, 1, 1.5, 2, 3, 4, 5, 6, 8, 12. To present the whole set of the boundaries 
on the same figure, each boundary is rescaled by a corresponding factor 
exp(— r/2). We would like to note here that the effect of rescaling is equiv- 
alent to the non- rescaled process with drift n = 1. One can see how the 
original ellipse-like shape (innermost figure) evolves through a twisted-like 
shape to the set of double-ring structures. The inner ring, always containing 
the origin, is so small on the scale of the Fig. [7| that is not visible. At r = 4, 
corresponding to the curve with an inner loop, we observe a topological phase 
transition. The support of the spectrum is no longer simply connected, for 
r > 4 it is annulus-like, i.e. eigenvalues are expelled from the central region. 
For even larger times, the outside rim of the annulus approaches the circle, 
and the inner one shrinks to the point z = in the t = oo limit. Indeed for 
large r the radius of the inner ring behaves like 



The outer boundary then forms approximately a circle with the radius e"^/^. 
To visualize the repulsion of the eigenvalues from the region around z = we 
performed higher statistics simulations for r = 4 and r = 5, corresponding 
to Figs I9|1(H respectively. Note that the eigenvalue distributions are quite 
high for small z and the size of the inner ring is indeed very small. Therefore 
it is quite difficult to observe numerically the exact exclusion of eigenvalues 
from the marked circle. 

Figure |H1 shows the comparison of numerically generated spectra versus 
the analytical prediction of the shape of the support of eigenvalues of the 



8 Numerical Tests 




(100) 
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Figure 7: Evolution of the rescaled (see text) contour of non-holomorphic do- 
main in the eigenvalue spectrum of the Ginibre-Girko multiplicative diffusion 
function of several times r. 

"evolution operator". Again, the same rescaling by exp(— r/2) was applied. 

Figure ^2 shows the evolution of the boundary for several sample evo- 
lution times r = 0.5,1,1.5,2,3,4,5,6,8,12 in the case of GUE. Again, to 
present the whole set of boundaries on the same figure, each boundary was 
rescaled by a corresponding factor exp(— r/2). One can observe how the 
almost real (for small times r) long-cigar shaped spectrum evolves into a 
broader shape, developing again singular behavior at r = 4 at the leftmost- 
edge of the spectrum. Again, the spectrum develops a topological transition, 
although it is not as clearly visible in the figure as in the complex case. The 
reason is that inversion symmetry has an additional exponential suppression 
factor, i.e. r — > e~'^/r. For very large times, the outside rim of the spectrum 
approaches the circle, and the inner one shrinks to a point in the t —>■ oo 
limit. 

As in the complex case, here we also show a sample simulation (Fig. I12|l 
versus the analytical prediction for the shape of the island. 

Finally in figures El and IHl above, we demonstrated the "localization- 
delocalization" phase transition in the case of the product of two hermitian 
matrix models of the considered type. Below Tcrit, a priori complex spectra 
condense on the positive real axes. At r = 1/2, the eigenvalues start to 
diffuse onto the x < half-plane. For very large times r, eigenvalues start 
to appear on both sides of the x = axis. 
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Figure 8: Comparison of the analytical contour for the eigenvalue spectrum 
of Ginibre-Girko multiplicative diffusion at evolution time r = 1., versus the 
numerical simulation of the spectrum. The generated ensemble consisted of 
60 matrices, each for N = M = 100. Note that the vertical axis is located at 
X = 1 and not at x = 0. The origin lies outside the figure. 




Figure 9: Comparison of the analytical contour for the eigenvalue spectrum 
of Ginibre-Girko multiplicative diffusion at evolution time r = 4., versus the 
high-statistics numerical simulation of the spectrum. 
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Figure 10: Comparison of the analytical contour for the eigenvalue spectrum 
of Ginibre-Girko multiplicative diffusion at evolution time r = 5, versus the 
numerical simulation of the spectrum. The presence of few eigenvalues inside 
the inner contour is caused by numerics, i.e. finite N and finite M effects. 



Figure 11: Evolution of the contour of non-holomorphic domain in the eigen- 
value spectrum of the GUE multiplicative diffusion as a function of evolution 
times r = 0.5, 1, 1.5, 2, 3, 4, 5, 6, 8, 12, from inside to outside, respectively. 
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Figure 12: Comparison of the analytical contour for the the eigenvalue spec- 
trum of GUE multiplicative diffusion at evolution time r = 1., versus the 
numerical simulation of the spectrum. The ensemble consists of 60 matrices, 
each obtained as a string of M = 100 of by matrices, where N = 100. 



9 Stability and Lyapunov Exponents 

In this section, we briefly discuss the relation between the spectrum of the 
"evolution operator" Ym and Lyapunov spectra. In the studies of products 
of random matrices one of the central issues is the use of limit theorems 
that would provide an information about the asymptotic limit of the rates 
of growth and of the spectrum of the product for a sequence of matrices. In 
particular the limit 

hm l,(ln|yM|) = Ai (101) 

defines the maximum Lyapunov characteristic exponent. The existence of 
this limit is assured for an infinite product of matrices by the Furstenberg [T^ 
theorem 

M 

YM = Y[^i^ (102) 

1=1 

for independent random matrices X{i) characterized by the probability dis- 
tribution function (i/i[X] = 'P[X](i[X]. It can be shown under very general 
assumptions that Ai is a non-random, self- averaging quantity, so that the 
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brackets (. . .) in the above expression can be dropped in almost any realiza- 
tion. The best known example is the use of Ai in the description of the DC 
conductivity of a one dimensional disordered system coupled to two electron 
reservoirs. The conductivity /t of the system, given by the Landauer formula 
yields for large M 

K ^ (103) 

where Ai is the maximal Lyapunov exponent of the product Ym- By the Bor- 
land conjecture, the inverse of Ai gives the localization length measuring the 
decay of an eigenvector of the system's Hamiltonian over the chain composed 
of M units. 

It turns out that the formalism of generalized Green functions (Section 
2 and 3) is particularly suitable for deriving an analytical formula EH] 
for for the one-dimensional disordered wires within the tight-binding 
Hamiltonian with a site diagonal disorder. Extensions to generalized Lya- 
punov exponents L{q) [13^ related to the exponential growth rates of the 
moments of the matrix product 

L(g)= hm ^Indn/I^) (104) 

have been also reported for systems described by Markovian evolutionary 
operators f36] used in the study of DC conductance statistics of the ran- 
dom dimer model introduced for explaining the exceptionally high electronic 
transport properties of some conjugated polymer chains. 

The eigenvectors of the product V^/^a/ form an orthonormal set of vectors 
(Lyapunov basis) corresponding to characteristic Lyapunov exponents whose 
exponentials, exp(Aj) are eigenvalues of the matrix \yIiYm\^^'^^^ in the limit 
of large M (Oseledec's theorem). 

The Lyapunov exponents are closely related to the linear stability analysis 
of the dynamic system x = F{x) for which small perturbations y around a 
given solution x{t) evolve in time according to ^(^2) = K{t2,ti)y{ti). They 
are defined as 

A, = hm ^^log\\K{t2,h)Mh)\\ (105) 

and for ergodic systems their spectrum describes global properties of the 
system's attractor (does not depend on the initial point). In contrast to the 
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Lyapunov spectra, the Lyapunov basis for the ergodic sequence of matrices 
describes local properties of the PRM system. 

Here we would like to underline, that Lyapunov exponents are in general 
different from the spectrum of eigenvalues obtained in this paper. Since the 
first Lyapunov exponent is defined as a property of the norm of the operator 
y, it is related to the spectral properties of the Y'^Y operator. This is why 
the Lyapunov exponents are always real. The spectra studied in this paper 
correspond to so-called stability exponents, which are directly related to the 
eigenvalues of the Y operator itself. As we have demonstrated above, the 
eigenvalues of Ym are in general complex numbers and so are the stability 
exponents: 

a,= hm ^(InA,) (106) 

Their real parts coincide with the Lyapunov exponents Aj = Re{ai) only if 
the matrix Ym can be written in a diagonal form. 

The relation between the stability and Lyapunov exponents in a broad 
class of dynamical systems has been investigated by Goldhirsch et al. ^7] . 
Motivated by mostly numerical results, the authors have conjectured that 
Lyapunov exponents equal stability ones for the ergodic dynamical systems 
whose spectrum is non degenerate and the attractor is bounded. The geo- 
metrical argument behind that statement is that generically, stretching and 
shrinking of the relative separation of "trajectories", starting at a given ini- 
tial point happen along fixed directions in the phase space. Taken from 
that perspective further elucidation would be required to clarify the relation 
between the Lyapunov stability and spectral properties of the generalized 
stability matrix investigated in this work. 



10 Summary 

We have introduced a natural generalization of the concept of geometric ran- 
dom walk ('geometric diffusion') in the space of large, non-commuting matri- 
ces. In particular, we proposed a simple diagrammatic method, allowing to 
calculate spectral properties of the resulting evolution operators as functions 
of the "evolution" time r. The construction was explicitly presented in the 
case of complex Gaussian ensembles and CUE ensembles. 

Our method was based on an exact relation between the eigenvalues of a 
product matrix and the eigenvalues of a certain block matrix. The essential 
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simplification is the fact tliat tlie random matrices appear linearly in tlie 
block matrix. Although we used the above identification only in the large 
limit, it could a-priori be used to study the properties of the spectrum of the 
product matrix also in the microscopic limit, i.e. on the scale of eigenvalue 
spacing. But then different, non-diagrammatic, methods would have to be 
employed. 

Despite the simplicity of the random matrix models chosen as the build- 
ing blocks of the matrix-valued diffusion, the spectrum develops a surprising 
feature, namely a region without eigenvalues appears within the spectrum, 
thus changing its topological properties. We can thus describe this behavior 
as a topological phase transition. This points at the appearance of partic- 
ular repulsion mechanism for sufficiently large evolution times. Our results 
were based on the simulation based on Gaussian ensembles. We are how- 
ever tempted to conjecture, that unfolding the spectrum in the vicinity of 
the point where the topological change occurs may unravel a possible novel 
universal spectral behavior. We would like to mention, that the critical be- 
havior described here as a "topological phase transition" appeared only in 
the limit M — > cxd, as a consequence of the inversion symmetry, acquired in 
this limit, of the equation describing the boundary of the two-dimensional 
support of eigenvalues. We would like to mention, that similar, qualitative 
effect of expulsion of complex eigenvalues from the certain region of com- 
plex plane was observed by Feinberg and Zee (first reference in j2H])- They 
considered however, a single random complex ensemble with non-Gaussian 
measure and they imposed the rotational symmetry of the spectrum. In 
our case, we considered the infinite product representing the Ginibre-Girko 
(therefore Gaussian) multiplicative diffusions, and we obtained the structural 
change of the spectrum for the general case, with non-trivial dependence on 
the polar variables. It is an intriguing problem for further studies, to what 
extent these results reflect the general "single ring hypothesis". 

We pointed out, in agreement with |37j, that stability exponents consid- 
ered in our paper carry much richer information about the properties of the 
system than the Lyapunov exponents. 

As a by-product of our analysis, we found a previously unknown (as far 
as we know) critical behavior in the case of the particular product of two 
hermitian ensembles. The a-priori complex spectra of the product are nev- 
ertheless real for evolution times below a certain critical value Tcrit, and 
migrate into the complex plane only for r > Tcru- A similar in spirit 
localization-delocalization phase transition was recently observed for a class 
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of non-hermitian Anderson Hamiltonians |38j . 

One of the motivations for this work was to propose a general formahsm, 
which can provide a straightforward method of analyzing spectral properties 
of multivariate diffusion- like processes, with the idea that our method could 
be used in different branches of theoretical physics and interdisciplinary ap- 
plications. Some particular applications will be presented elsewhere. Finally, 
we would like to mention, that the presented formalism allows also to estab- 
lish a direct link to the diffusion processes based on Free Random Variables 
techniques 1^ 1^ , in particular it can demonstrate the emergence of the 
complex Burgers equations governing the evolution of the spectral generating 
functions 
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